#ifndef RSIM_EIG_H
#define RSIM_EIG_H

/* -------------------------------------------------------------------------- *
 * File: Eig.h							      	      *
 * Authors: Kishor Bhalerao                                                   *
 * Email : kishor8dm@gmail.com                                                *
 * Contributors:                      					      *
 * Email: 								      *
 * copyright (c) 2010 Authors.                                                *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify it    *
 * under the terms of the GNU General Public License as published by the Free *
 * Software Foundation, either version 3 of the License, or any later version.*
 *									      *
 * This program is distributed in the hope that it will be useful, but        *
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
 * for more details. (http://www.gnu.org/licenses/)			      *
 * -------------------------------------------------------------------------- */

#include "LinAlgDecl.h"
#include "Matrix.h"
#include "Vector.h"

#include <algorithm>
#include <cmath>

using namespace std;

namespace RSIM
{

/** 
	
    Note: This file was obtained from Template Numerical Toolkit (TNT).
    Details of TNT can be found on the website http://math.nist.gov/tnt/
    This file is minimally modified to make it compatible with RSIM.
    The originial file name (available for download on TNT website) 
    is jama_eig.h.

<P>
    Computes eigenvalues and eigenvectors of a real (non-complex)
    matrix. 
<P>
    If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
    diagonal and the eigenvector matrix V is orthogonal. That is,
	the diagonal values of D are the eigenvalues, and
    V*V' = I, where I is the identity matrix.  The columns of V 
    represent the eigenvectors in the sense that A*V = V*D.
    
<P>
    If A is not symmetric, then the eigenvalue matrix D is block diagonal
    with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
    a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
    eigenvalues look like
<pre>

          u + iv     .        .          .      .    .
            .      u - iv     .          .      .    .
            .        .      a + ib       .      .    .
            .        .        .        a - ib   .    .
            .        .        .          .      x    .
            .        .        .          .      .    y
</pre>
        then D looks like
<pre>

            u        v        .          .      .    .
           -v        u        .          .      .    . 
            .        .        a          b      .    .
            .        .       -b          a      .    .
            .        .        .          .      x    .
            .        .        .          .      .    y
</pre>
    This keeps V a real matrix in both symmetric and non-symmetric
    cases, and A*V = V*D.
    
    
    
    <p>
    The matrix V may be badly
    conditioned, or even singular, so the validity of the equation
    A = V*D*inverse(V) depends upon the condition number of V.

   <p>
	(Adapted from JAMA, a Java Matrix Library, developed by jointly 
	by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).
**/

class Eigenvalue{
	private:
		
		/** Row and column dimension (square matrix).  */
		int n;

		int issymmetric; /* boolean*/

		/** Arrays for internal storage of eigenvalues. */

		Vector d;         /* real part */
		Vector e;         /* img part */

		/** Array for internal storage of eigenvectors. */
		Matrix V;

		/** Array for internal storage of nonsymmetric Hessenberg form.
		@serial internal storage of nonsymmetric Hessenberg form.
		*/
		Matrix H;
		
		/** Working storage for nonsymmetric algorithm.
		@serial working storage for nonsymmetric algorithm.
		*/
		Vector ort;

		// Symmetric Householder reduction to tridiagonal form.
		void tred2();

		// Symmetric tridiagonal QL algorithm.
		void tql2 ();

		// Nonsymmetric reduction to Hessenberg form.
		void orthes ();

		// Complex scalar division.
		double cdivr, cdivi;
		void cdiv(double xr,double xi,double yr,double yi);

		// Nonsymmetric reduction from Hessenberg to real Schur form.
		void hqr2();
	public:

		/** Check for symmetry, then construct the eigenvalue decomposition
		@param A    Square real (non-complex) matrix */
		Eigenvalue(const Matrix& A);

		/** Return the eigenvector matrix
		@return    V */
		void getV(Matrix& V_){
			V_ = V;
			return;
		}

		/** Return the real parts of the eigenvalues
		@return     real(diag(D)) */
		void getRealEigenvalues(Vector& d_){
			d_ = d;
			return ;
		}

		/** Return the imaginary parts of the eigenvalues
		in parameter e_.
		@param e_: new matrix with imaginary parts of the eigenvalues. */
		void getImagEigenvalues (Vector& e_) {
			e_ = e;
			return;
		}

   
/** 
	Computes the block diagonal eigenvalue matrix.
    If the original matrix A is not symmetric, then the eigenvalue 
	matrix D is block diagonal with the real eigenvalues in 1-by-1 
	blocks and any complex eigenvalues,
    a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
    eigenvalues look like
<pre>

          u + iv     .        .          .      .    .
            .      u - iv     .          .      .    .
            .        .      a + ib       .      .    .
            .        .        .        a - ib   .    .
            .        .        .          .      x    .
            .        .        .          .      .    y
</pre>
        then D looks like
<pre>

            u        v        .          .      .    .
           -v        u        .          .      .    . 
            .        .        a          b      .    .
            .        .       -b          a      .    .
            .        .        .          .      x    .
            .        .        .          .      .    y
</pre>
    This keeps V a real matrix in both symmetric and non-symmetric
    cases, and A*V = V*D.

	@param D: upon return, the matrix is filled with the block diagonal 
	eigenvalue matrix.
	
*/
		void getD (Matrix& D);
}; // class Eigenvalue
} //namespace RSIM

#endif
// RSIM_EIG_H
